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Abstract 

In  this  paper,  we  examine  alternative  methods  for  re¬ 
ducing  the  dimensionality  of  nonlinear  dynamical  sys¬ 
tem  models  arising  in  control  of  rapid  thermal  chemical 
vapor  deposition  (RTCVD)  for  semiconductor  manufac¬ 
turing.  We  focus  on  model  reduction  for  the  ordinary 
differential  equation  model  describing  heat  transfer  to, 
from,  and  within  a  semiconductor  wafer  in  the  RTCVD 
chamber.  Two  model  reduction  approaches  are  stud¬ 
ied  and  compared:  the  proper  orthogonal  decomposition 
and  the  method  of  balancing.  This  leads  to  a  discussion 
of  computational  issues  in  the  practical  implementation 
of  balancing  for  nonlinear  systems. 

1  Introduction 

Model  reduction  deals  with  methods  for  reducing  the 
dimensionality  of  dynamical  system  models.  The  moti¬ 
vation  is  that  models  of  lower  dimension  are  less  com¬ 
plex  and  easier  to  work  with  for  various  purposes  such 
as  simulation,  optimization,  and  control.  One  source  of 
dynamical  system  models  which  are  good  candidates  for 
model  reduction  is  rapid  thermal  chemical  vapor  depo¬ 
sition  (RTCVD),  a  process  used  to  deposit  thin  films  on 
a  semiconductor  wafer  via  thermally  activated  chemical 
mechanisms. 

This  paper  describes  our  recent  progress  and  on¬ 
going  investigation  in  nonlinear  model  reduction  for 
RTCVD.  These  efforts  have  taken  place  as  part  of  a  joint 
industry-academia  research  project  to  optimize  the  epi¬ 
taxial  growth  of  silicon-germanium  (Si-Ge)  heterostruc¬ 
tures  in  a  commercial  RTCVD  reactor.  Participants  are 
the  University  of  Maryland’s  Institute  for  Systems  Re¬ 
search  (ISR)  of  College  Park,  MD,  and  the  Northrop 
Grumman  Corporation’s  Electronic  Sensors  and  Systems 
Division  (ESSD)  of  Linthicum,  MD.  Details  and  some 
results  from  this  project  are  described  in  [1,  2]. 

Our  approach  to  model  development  for  control  and 
optimization  combines  first-principles,  experimental  and 
simulation  data,  and  system-theoretic  ideas  for  purposes 
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of  model  validation  and  model  reduction.  This  paper 
focuses  on  two  main  topics:  a  comparison  of  model 
reduction  approaches  as  applied  to  one  component  of 
the  overall  RTCVD  model,  and  computational  issues 
in  nonlinear  balancing.  Attractive  properties,  draw¬ 
backs,  and  challenges  of  model  reduction  via  established 
statistically-based  empirical  methods  and  balanced  real¬ 
ization  methods  are  discussed,  demonstrated,  and  com¬ 
pared  via  numerical  studies. 

2  RTCVD  Models 

Process  and  equipment  models  for  RTCVD  involve  pre¬ 
diction  of  time  evolution  for  heat  transfer  within  and 
among  the  wafer,  chamber  walls,  and  process  gases  (tem¬ 
perature  fields  in  solids  and  gases),  momentum  transport 
in  the  gas  phase  (gas  flow  velocity  vectors),  mass  trans¬ 
port  in  the  gas  phase  (species  concentrations) ,  and  chem¬ 
ical  reaction  kinetics  in  the  gas  phase  and  at  the  wafer 
surface  (reaction  rates,  deposition  thickness).  Thus, 
they  consist  mainly  of  balance  equations  for  conservation 
of  energy  (heat  equation),  momentum  (Navier-Stokes), 
and  mass  (continuity),  along  with  equations  that  de¬ 
scribe  the  relevant  chemical  mechanisms  (e.g.  Arrhenius 
laws)  and  boundary  and  initial  conditions. 

In  their  continuum  form,  the  balance  equations  for 
RTCVD  yield  a  system  of  nonlinear  coupled  partial  dif¬ 
ferential  equations  (PDEs)  with  associated  boundary 
conditions  (BCs)  and  initial  conditions  (ICs).  Lumped 
versions  of  the  equations  can  be  obtained  using  an  ap¬ 
propriate  discretization  method,  e.g.,  finite-elements,  to 
yield  a  system  of  coupled  nonlinear  ordinary  differen¬ 
tial  equations  (ODEs).  The  system  of  ODEs  can  be 
decoupled  by  invoking  certain  simplifying  assumptions. 
Even  the  resulting  simplified  nonlinear  system  is  typ¬ 
ically  of  relatively  high-order,  so  that  not  only  is  the 
model  computationally  demanding  for  simulation,  but 
moreover  it  is  computationally  prohibitive  for  real-time 
control.  Thus,  the  motivation  for  reducing  the  model 
order  is  apparent. 

To  illustrate  ideas  in  model  reduction,  we  focus  on  the 
evolution  system  describing  heat  transfer  on  the  surface 
of  a  silicon  wafer  in  the  ASM  Epsilon-1,  a  commercial 
RTCVD  reactor  used  by  Northrop  Grumman  ESSD.  The 
Epsilon-1  has  a  horizontally  oriented  process  chamber, 
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where 

[F]ij  =  [Ar]ij  r| 

and  x  and  u  are  translations  of  Tw  and  P,  respectively. 


3  Model  Reduction 


Figure  1:  Semiconductor  wafer  discretized  into  annular 
regions.  Heat  transfer  in  wafer  is  described  by  energy 
balance  for  each  discrete  element. 

with  gases  flowing  parallel  to  and  over  the  top  surface  of 
the  wafer.  The  silicon  wafer  is  very  thin  and  rotating  on 
a  susceptor  so  that  it  is  reasonable  to  assume  axial  and 
azimuthal  symmetry,  i.e.,  all  mechanisms,  including  heat 
transfer,  are  functions  only  of  time  t  and  radial  position 
r  on  the  wafer  surface.  The  wafer  is  heated  by  lamps, 
grouped  into  “lamp  zones” ,  whose  radiant  intensity  pro¬ 
files  are  functions  of  r  with  corresponding  power  levels 
that  are  available  as  control  inputs.  The  relevant  heat 
transfer  mechanisms  are  shown  in  Figure  1. 

A  simplified  version  of  the  temperature  field  evolution 
on  the  surface  of  a  silicon  wafer  during  RTCVD  is  given 
by  the  ODE 

Tw  =  Ac  Tw  +  Ar  +  Av  Tw  +  T  +  B  P  (1 ) 

where  Tw(t)  is  a  n-vector  of  temperatures  correspond¬ 
ing  to  radial  positions  on  the  wafer  surface,  Ac,  Ar,  and 
Av  are  constant  nxn  matrices  representing  the  effects  of 
conductive,  radiative,  and  convective  heat  transfer  mech¬ 
anisms,  respectively,  F  is  a  constant  n-vector  which  ac¬ 
counts  for  the  gas  and  chamber  wall  temperature,  B  is  a 
n  x  m  matrix  of  discretized  lamp  zone  radiant  intensity 
profiles,  and  P  is  a  m-vector  of  control  inputs  corre¬ 
sponding  to  lamp  zone  power  levels.  Note  that  for  the 
nonlinear  term  the  exponent  is  taken  component¬ 
wise.  Our  model  uses  m  =  3  lamp  zone  actuators  and 
n  =  101  grid  points  for  the  spatial  discretization.  The 
parameters  and  lamp  zone  intensity  profiles  used  in  Ac, 
Ar.  Av,  F,  and  B  are  derived  from  material  properties, 
radiative  heat  transfer  analysis,  and  experimental  results 
(see  [2]). 

To  model  the  measurement  of  temperature  at  discrete 
points  on  the  wafer  surface  via  thermocouples,  we  aug¬ 
ment  the  nonlinear  state  equation  (1)  with  the  linear 
output  equation 

Ttc  =  CTw  (2) 

where  Ttc  is  a  p-vector  of  thermocouple  measurements, 
and  C  is  a  p  x  n  matrix  with  entries  corresponding  to 
thermocouple  locations.  We  use  p  =  3  thermocouples  in 
our  model,  at  the  center,  edge,  and  midpoint  between 
center  and  edge. 

Later,  we  make  use  of  a  linearized  version  of  (1).  The 
state  equation  is  linearized  about  the  constant  tempera¬ 
ture  profile  Tw  =  F  and  is  of  the  form 

(3) 


A  general  approach  to  model  reduction  is  to  find  a  coor¬ 
dinate  transformation  of  the  state  space  under  which  the 
state  components  can  be  ranked  in  a  meaningful  way  in 
terms  of  their  influence  on  the  system  behavior.  Then, 
state  components  of  the  transformed  system  with  rela¬ 
tively  small  influence  can  be  truncated  without  substan¬ 
tially  degrading  the  correctness,  i.e.,  predictive  capabil¬ 
ity,  of  the  model.  We  note  that  for  systems  evolving 
on  1R",  each  coordinate  transformation  can  be  identified 
with  a  corresponding  set  of  basis  n-vectors. 

Proper  Orthogonal  Decomposition 

One  approach  to  finding  a  basis  for  the  desired  co¬ 
ordinate  transformation  is  application  of  the  Karhunen- 
Loeve  decomposition  (see,  e.g.,  [3,  4,  5,  6]),  also  known  as 
the  proper  orthogonal  decomposition  (POD),  method  of 
empirical  eigenfunctions,  or  principal  components  anal¬ 
ysis  (PCA).  The  POD  is  a  statistical  pattern  analysis 
technique  for  finding  the  dominant  structures  in  an  en¬ 
semble  of  spatially  distributed  data.  These  structures 
can  be  used  as  an  orthogonal  basis  for  efficient  repre¬ 
sentation  of  the  ensemble.  In  particular,  the  POD  basis 
elements  are  the  orthogonal  eigenfunctions  of  the  two- 
point  spatial  covariance  of  the  data  ensemble.  When  the 
data  ensemble  consists  of  vectors  in  1R",  the  POD  basis 
vectors  are  just  the  the  columns  of  the  matrix  U  in  the 
singular  value  decomposition 

X  =  UEVt  (4) 

where  X  is  a  matrix  whose  columns  are  the  members  of 
the  data  ensemble. 

For  the  case  of  a  dynamical  system  describing  a  flow, 
e.g.,  a  heat  flow,  the  data  ensemble  is  typically  time  se¬ 
ries  data,  i.e.,  “snapshots”  of  the  flow  captured  at  equally 
spaced  intervals  in  time.  The  POD  coordinate  transfor¬ 
mation  diagonalizes,  or  decouples,  the  covariance  of  the 
time  series  data.  The  basis  elements  are  the  principal 
axes  of  the  flow  which  generated  the  time  series  empiri¬ 
cal  data.  Each  has  a  corresponding  eigenvalue  (given  by 
the  entries  on  the  diagonal  of  E  in  (4)),  which  provides  a 
measure  of  the  relative  “energy”,  i.e.,  mean  square  fluc¬ 
tuation,  associated  with  that  particular  direction  in  the 
state  space.  This  measure  can  also  be  interpreted  as  the 
relative  amount  of  time  that  the  flow  spends  along  the 
corresponding  principal  axis.  Thus,  it  serves  as  a  well- 
defined  measure  of  the  influence  of  a  state  component  on 
the  system  behavior. 

The  POD  basis  is  optimal  from  the  points  of  view  of 
data  compression  and  error  minimization.  Specifically,  a 
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Figure  2:  Lamp  power  settings  for  RSC  recipe. 


truncated  series  representation  of  the  data  has  a  smaller 
mean  square  error  than  a  representation  by  any  other 
basis  of  the  same  dimension  (see,  e.g.,  [3]). 


0.2 

0.1 

0 

-0.1 

-0.2 


0  0.2  0.4  0.6  0.8  1 

0.03% 


Figure  3:  Basis  elements 
RSC  empirical  data. 
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Application 

The  attractive  properties  of  the  POD  have  led  to  suc¬ 
cess  in  applying  it  to  areas  such  as  turbulence  model¬ 
ing  (e.g.,  [3])  and  pattern  recognition  (e.g.,  [7]).  Re¬ 
cently,  much  work  has  been  done  to  study  its  use  for 
RTCVD  model  reduction  (e.g.,  [8,  9,  10,  11]).  We  have 
performed  a  similar  study,  applying  the  POD  technique 
toward  finding  a  low-order  approximation  to  (1),  de¬ 
scribed  as  follows. 

To  generate  empirical  time  series  data,  i.e. ,  snapshots 
of  the  wafer  temperature  field,  the  system  (1-2)  was  sim¬ 
ulated  using  two  different  types  of  control  input  recipes. 
They  are  referred  to  as  Ramp-Soak-Cool  (RSC)  and 
Perturbation-Of-Constant  (POC). 

The  RSC  recipe  mimics  a  typical  processing  recipe  in 
which  a  lamp  zone  power  setting  is  gradually  ramped 
up  from  zero  to  full  power,  maintained  at  full  power  for 
a  specified  period  of  time,  and  then  gradually  ramped 
down  from  full  to  zero  power,  as  shown  in  Figure  2.  This 
recipe  is  applied  to  one  of  the  lamp  zones,  while  the 
other  two  zones  are  held  at  zero  power.  The  simulation 
is  then  repeated  using  RSC  individually  for  each  of  the 
other  two  lamp  zones.  The  entire  ensemble  (three  sets) 
of  time  series  data  is  combined  and  used  to  compute  the 
POD  basis  elements,  which  are  then  ranked  according  to 
magnitude  of  associated  eigenvalue.  The  basis  elements 
with  the  four  largest  eigenvalues  are  shown  in  Figure  3. 
Corresponding  relative  energy  values  are  contained  in 
Table  1. 

The  POC  recipe  applies  small  perturbations  of  a  pre¬ 
determined  set  of  constant  power  settings  which,  if  left 
unperturbed,  would  result  in  a  uniform  steady  state  tem¬ 
perature  field  of  1000K.  The  perturbations  are  achieved 
by  adjusting  the  power  setting  of  each  lamp  zone,  one 
at  a  time,  first  to  110%  and  then  to  90%,  of  the  original 
setting.  This  results  in  6  different  control  recipes  (see  [2] 
for  details). 

The  system  response  to  excitation  from  each  of  the  six 
POC  recipes  is  sampled  and  combined  as  the  time  series 
data  for  computing  POD  basis  elements.  Once  again, 


Original:  n  =  101  Reduced:  k=1 


Figure  4:  Thermocouple  readings  for  original  and  re¬ 
duced  models  with  Test  Recipe  2  using  POD  RSC  coor¬ 
dinate  transformation.  Original  n  =  100.  Reduced  k  = 
1,  2,  3. 


basis  elements  are  ranked  by  corresponding  relative  en¬ 
ergy  value,  contained  in  Table  1. 

We  now  test  the  efficiency  of  each  POD  basis  by 
projecting  (1)  onto  a  fc-dimensional  subspace  spanned 
by  the  k  most  important  basis  vectors  using  standard 
Galerkin  projection,  where  k  <C  n.  The  reduced  or¬ 
der  approximation  is  numerically  integrated  using  two 
test  recipes  as  control  inputs,  that  are  different  from  the 
recipes  used  to  generate  the  RSC  and  POC  data  ensem¬ 
bles. 

Test  Recipe  1  P  =  [0.5  0.5  0.5]  0  <  t  <  60 

Test  Recipe  2  P  =  [1.0  0.0  0.0]  0  <  t  <  24 

P  =  [0.0  1.0  0.0]  24  <  t  <  42 

P  =  [0.0  0.0  1.0]  42  <  t  <  60 

Figure  4  shows  simulated  thermocouple  readings  for 
Test  Recipe  2,  where  the  k  basis  vectors  for  the  Galerkin 
projection  are  from  POD  RSC,  for  k  =  1,  2,  3.  For 
complete  results  see  [2],  The  error  between  original  and 
reduced  models  is  computed  as  the  maximum  deviation 
between  actual  and  approximate  thermocouple  readings. 
Results  for  all  test  simulations  are  contained  in  Table  2. 


Shortcomings 

It  is  clear  that  the  efficiency  of  a  basis  determined  via 
the  POD  method  depends  strongly  on  how  well  the  data 
ensemble  captures  the  relevant  system  behavior.  This 
leads  to  serious  shortcomings  for  model  reduction  of  con¬ 
trol  systems  with  inputs  and  outputs.  The  POD  basis 
elements  will  be  sensitive  to  the  choice  of  control  in¬ 
puts,  ICs,  and  BCs  used  to  generate  the  empirical  time- 
series  data.  For  nonlinear  systems,  small  perturbations 
in  these  parameters  may  produce  qualitatively  different 
system  responses.  In  addition,  the  data  may  fail  to  cap¬ 
ture  dynamical  effects  occuring  at  widely  differing  time 
scales.  Usually,  these  issues  are  ignored,  because  the 
model  is  only  being  used  for  a  particular  purpose,  e.g., 
to  simulate  tracking  of  one  particular  reference  trajec¬ 
tory.  In  that  case,  a  representative  set  of  control  inputs, 
BCs,  and  ICs  is  selected  for  generating  the  data  ensem¬ 
ble.  But  the  optimality  property  of  the  POD  basis,  and 
the  predicitive  capability  of  the  resulting  reduced  order 
model,  may  be  localized  to  a  relatively  small  region  of 
the  space  of  allowable  inputs,  BCs,  and  ICs.  In  addi¬ 
tion,  the  POD  approach  fails  to  consider  the  influence  of 
state  components  on  the  system  measurements,  or  out¬ 
puts.  It  would  be  desirable  to  truncate  state  compo¬ 
nents  that  have  small  influence  on  the  outputs,  i.e.  that 
do  not  appear  in  our  measurements.  The  POD  method 
does  not  do  this,  thus  its  efficiency  may  be  diminished 
when  constructing  “black-box”  input-output  models  of 
the  system. 

Balancing 

In  response  to  these  and  other  concerns,  we  considered 
an  alternative  approach  based  on  the  method  of  balanced 
realizations  (see,  e.g.,  [12,  13]).  In  this  method,  a  coor¬ 
dinate  transformation  is  computed  which  allows  state 
components  to  be  ranked  according  to  their  influence  on 
the  input-output  behavior  of  the  system  as  measured 
by  the  Hankel  norm  of  the  system,  i.e.,  the  gain  from 
past  inputs  to  future  outputs.  The  resulting  basis  for 
the  linear  transformation  makes  the  transformed  realiza¬ 
tion  “equally  controllable  and  observable,”  i.e.  balanced, 
and  depends  only  upon  intrinsic  properties  of  the  orig¬ 
inal  model,  specifically  controllability  and  observability, 
embodied  in  its  evolution  equation  and  output  equation. 
In  the  linear  case,  explicit  bounds  can  be  computed  for 
the  error  between  the  original  and  reduced  order  mod¬ 
els.  These  error  bounds  are  independent  of  any  partic¬ 
ular  set  of  control  input  signals,  BCs,  or  ICs.  Although 
explicit  error  bounds  have  not  yet  been  found  for  the 
nonlinear  case,  we  still  wish  to  exploit  the  property  that 
the  correctness  of  an  approximation  using  a  truncated 
balanced  realization  does  not  depend  upon  generating 
a  representative  data  ensemble.  Furthermore,  the  bal¬ 
ancing  method  emphasizes  state  components  that  are 
both  strongly  controllable  and  strongly  observable,  so 
that  state  components  which  are  least  likely  to  influence 
the  measurements  are  truncated. 


Application 

We  have  applied  the  balancing  method  to  the  lin¬ 
earized  control  system  given  by  (3)  and  (2).  The  lin¬ 
earized  version  is  used  because  efforts  toward  develop¬ 
ment  of  useful  algorithms  for  nonlinear  balancing  are 
still  underway,  as  described  later  in  this  paper. 

Numerical  difficulties  arise  when  applying  standard 
linear  balancing  algorithms  to  our  linearized  system.  In 
particular,  the  system  is  nearly  nonminimal,  i.e.,  the 
condition  numbers  of  the  controllability  and  observabil¬ 
ity  Gramian  matrices  are  very  large  (see  [2]  for  details). 
Safanov  and  Chiang  [14]  have  proposed  a  method  for 
overcoming  such  difficulties,  based  on  the  ordered  Schur 
decomposition  of  the  product  of  the  controllability  and 
observability  Gramian  matrices.  We  have  used  their 
method  to  produce  coordinate  transformation  basis  vec¬ 
tors  and  fcth-order  reduced  models  for  the  linearized  sys¬ 
tem.  The  results  are  shown  in  Tables  1  and  2. 

Discussion 

Due  to  the  shape  of  the  lamp  zone  heat  flux  intensity 
profiles  and  the  smoothing  effect  of  the  diffusion  oper¬ 
ator,  the  evolution  of  the  wafer  temperature  field  does 
not  produce  especially  interesting  behavior,  e.g.,  spatial 
profiles  whose  fluctuations  from  the  mean  vary  substan¬ 
tially  in  the  mean  square  sense  from  the  initial  profile, 
assuming  the  initial  profile  is  relatively  smooth.  Thus, 
we  expect  little  difficulty  in  capturing  the  essence  of  the 
input-output  behavior  of  the  system  in  a  low  dimen¬ 
sional  model.  Our  results  show  that  this  is  indeed  the 
case. 

For  the  inputs  used  in  the  validation  tests,  the  input- 
output  behavior  of  the  wafer  heat  transfer  system  can 
be  reconstructed  using  reduced  models  of  order  4  so 
that  thermocouple  readings  are  within  1  degree  C  of  the 
readings  obtained  using  the  original  model.  This  holds 
whether  the  POD  or  balancing  method  is  used,  and  for 
whichever  empirical  data  ensemble  was  used  for  comput¬ 
ing  the  POD  transformation.  Even  reduced  models  of  or¬ 
der  2  produce  a  reasonable  approximation  with  “worst 
case”  errors  less  than  15  degrees  C. 

The  POD  method  appears  to  have  performed  slightly 
better  than  the  balancing  method  in  this  study.  One 
reason  for  this  result  is  that  the  balancing  transforma¬ 
tion  was  computed  for  the  linearized  system,  while  the 
validation  tests  were  performed  for  the  reduced  order 
nonlinear  system.  Another  reason  is  the  simple  input- 
output  behavior  of  this  particular  system.  The  princi¬ 
pal  components  of  the  flow  are  relatively  insensitive  to 
the  choice  of  inputs,  and  hence,  any  set  of  empirically 
determined  eigenfunctions  for  this  system  will  likely  be 
efficient  for  model  reduction  purposes. 

Thus,  the  results  of  this  preliminary  study  are  not  de¬ 
cisive  regarding  choice  of  reduction  method.  We  seek  a 
model  reduction  approach  which  provides  explicit  error 
bounds,  does  not  depend  on  a  specific  data  ensemble, 
and  which  can  be  applied  directly  to  the  nonlinear  con¬ 
trol  system. 


Method 

Mode  1 

Mode  2 

Mode  3 

Mode  4 

POD  RSC 

95.06 

4.77 

0.14 

0.03 

POD  POC 

93.43 

6.25 

0.23 

0.09 

Balancing 

98.02 

1.83 

0.13 

0.02 

Table  1:  Normalized  eigenvalues,  i.e. ,  percent  energy, 
corresponding  to  basis  elements  used  in  model  reduction 
for  POD  method  with  RSC  data,  POD  method  with 
POC  data,  and  balancing  approach. 


Reduction 

Reduced  Model  Order 

Simulation 

Method 

1 

2 

3 

4 

Test  1 

POD  RSC 

27.23 

2.68 

0.58 

0.11 

POD  POC 

26.85 

1.26 

1.13 

0.10 

Balancing 

50.68 

7.03 

0.44 

0.08 

Test  2 

POD  RSC 

72.33 

5.22 

1.48 

0.18 

POD  POC 

72.60 

4.79 

4.35 

0.43 

Balancing 

80.81 

14.28 

1.70 

0.12 

Table  2:  Maximum  deviation  (degrees  C)  between  out¬ 
puts  of  original  and  reduced  models  for  POD  method 
with  RSC  data,  POD  method  with  POC  data,  and  bal¬ 
ancing  approach. 

4  Computational  Issues  In  Non¬ 
linear  Balancing 

In  order  to  address  some  of  the  deficiencies  in  model  re¬ 
duction  via  the  POD  and  linear  balancing  approaches, 
we  seek  a  procedure  to  apply  the  balancing  method  di¬ 
rectly  to  the  nonlinear  system  of  interest.  We  have  al¬ 
ready  seen  that  for  linear  systems,  the  balancing  coor¬ 
dinate  transformation  can  be  determined  efficiently  us¬ 
ing  known  algorithms.  The  theory  of  balancing  was  ex¬ 
tended  to  a  class  of  nonlinear  systems  by  Scherpen  [15], 
in  which  a  general  theory  and  procedure  is  proposed. 
In  contrast  to  the  linear  case,  the  theory  and  procedure 
for  nonlinear  balancing  is  not  immediately  amenable  to 
computational  implementation.  Here,  we  begin  to  ad¬ 
dress  some  important  computational  issues  in  balancing 
for  nonlinear  systems. 

We  consider  the  balancing  procedure  presented  in  [15] 
for  open  loop  stable  nonlinear  systems  of  the  form 

x  =  f(x)  +  g(x)  u  (5) 

V  =  h(x)  (6) 

where  x  =  (#i, . . . ,  xn)  £  Rn  are  local  coordinates  for  a 
smooth  state  manifold  denoted  by  M ,  u  £  !Rm,  y  £  1RP, 
/  and  h  are  of  class  C°°,  /( 0)  =  0,  and  h(0)  =  0. 

To  determine  the  balancing  state  space  transformation 
we  seek  a  change  of  coordinates  under  which  the  trans¬ 
formed  system  realization  is  equally  controllable  and  ob¬ 
servable.  The  degree  to  which  a  system  realization  is 
controllable  or  observable,  respectively,  can  be  measured 


in  a  precise  way  using  the  controllability  and  observabil¬ 
ity  energy  functions  of  the  system  [15]. 

Definition  4.1  The  controllability  function, 

Lc  :  1R”  — y  1R,  and  observability  function,  La  :  lRra  — >  It, 
for  (5-6)  are  defined  by 

1  hQ 

Lc{x 0)  =  min  -  /  \\u(t)\\2  dt  (7) 

U  £  £ 2(  — OO,0)  ^  j- oo 

x(—oo)  =  0 
*(0)  =  xo 

and 

1  r°° 

Lo(x0)  =  -J0  \\y{t)\\2dt,  (8) 

x(0)  =  xq  ,  u(t)  =  0 ,  0  <  t  <  oo. 

In  the  linear  case,  the  controllability  and  observability 
functions  specialize  to  the  quadratic  forms 

Lc(x 0)  =  ^Xq  WcT1Xo  (9) 

L0{x  0)  =  ^XqWoXo  (10) 

where  the  constant  n  x  n  matrices  Wc  and  W0  are  the 
familiar  controllability  and  observability  Gramian  ma¬ 
trices,  respectively. 

It  is  not  surprising,  then,  that  a  key  step  in  the  nonlin¬ 
ear  balancing  procedure  of  [15]  is  to  find  a  change  of  co¬ 
ordinates  under  which  Lc  is  locally  quadratic  in  a  neigh¬ 
borhood  U  of  its  critical  point  0.  This  is  accomplished 
by  appealing  to  the  Morse  Lemma  (see,  e.g.,  [16]). 

Theorem  4.2  (Morse  Lemma)  Let  p  be  a  non¬ 
degenerate  critical  point  for  the  real-valued  function  F 
with  index  r.  Then  in  a  neighborhood  U  of  p  there  exists 
a  local  change  of  coordinates  zi->i  =  'P(z),  'L(O)  =  p, 
such  that 

r  n 

F(x)  =  F(*{z))  =  F(p)  E  zi  (n) 

i—1  i=r+ 1 

holds  for  all  x  £U. 

There  are  various  difficulties  in  the  computation  of 
Lc  and  L0 ,  and  numerical  implementation  of  the  Morse 
Lemma.  We  are  currently  investigating  several  ap¬ 
proaches  for  each  problem.  Here,  we  discuss  the  prob¬ 
lem  of  diagonalizing  a  particular  ^-dependent  symmetric 
matrix  on  a  neighborhood  of  0,  which  arises  in  determi¬ 
nation  of  a  smooth  Morse  coordinate  transformation  for 
a  function  F  around  a  critical  point  0.  The  key  is  that 
the  diagonalization  must  be  done  in  a  smooth  way  on 
1R”,  and  thus  results  in  a  search  over  smooth  functions 
taking  values  in  the  matrix  Lie  group  SO(n). 

It  is  known  (see,  e.g.,  [16])  that  for  F  with  critical 
point  0  there  exists  a  symmetric  matrix  H( x)  for  each  x 
in  a  neighborhood  U  of  0  such  that 

F{x)  =  xT  Fl(x)  x 


(12) 


with  H( 0)  =  dfF{ 0).  Suppose  we  have  an  algorithm 
that  produces  the  desired  H{x)  for  each  x.  The  next 
step  in  implementation  of  the  Morse  Lemma  is  to  find  a 
smooth  diagonalization  of  the  matrix  H(x),  i.e. ,  find  a 
smooth  function  U(x)  taking  values  in  SO(n)  such  that 


UT  (x)  H  (x)  U  (x)  =  A(x)  =diag{Ai(a:),...,An(a:)} 

(13) 

for  each  x  in  U. 

We  propose  the  following  approach.  Choose  a  basis 
{ai,a2,  ■  ■ .}  for  the  £2-completion  of  CfffU,^),  i.e., 
smooth  real-valued  functions  with  compact  support  on 
the  closure  of  U ,  and  a  basis  {A\,...,Ak}  for  so{n), 
i.e.,  skew  matrices,  with  k  =  n(n  —  l)/2.  Then  a  Nk- 
dimensional  approximation  to  the  desired  U(x)  is  given 

by 

k  N  \ 

A‘  J2  cnaj(xn  (14) 

*=  1  3= i  / 


U  (x)  =  exp 


and  can  be  found  by  solving  the  finite-dimensional  min¬ 
imization  problem 


min  f  \\U(x)H(x)U(x)  —  A(x)\\dx.  (15) 

Cij  Ju 

i  =  1, . . . ,  k 
j  =  l,...,N 


5  Conclusion 

We  have  presented  a  comparison  of  nonlinear  model  re¬ 
duction  methods,  POD  and  balancing,  which  appear  to 
perform  well  for  a  specialized  problem  in  RTCVD,  but 
give  no  guarantees  of  satisfactory  performance  in  general 
for  dimensionality  reduction  of  nonlinear  systems  with 
control  inputs  and  outputs.  We  address  these  deficien¬ 
cies  by  appealing  to  the  theory  of  nonlinear  balancing, 
and  point  out  difficulties  in  terms  of  numerical  imple¬ 
mentation.  We  focus  on  one  step  involved  in  computa¬ 
tion  of  the  Morse  coordinate  transformation,  in  which 
a  smooth  5'0(?r)-valued  function  must  be  determined 
to  diagonalize  a  particular  x-dependent  symmetric  ma¬ 
trix.  We  suggest  an  approach  which  requires  solution  of 
a  finite-dimensional  minimization  problem. 
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